import numpy as np
import mpmath as mp

class sate_motion(object):
    def __init__(self, microthrust, sate_mass) -> None:
        '''
        计算卫星运动所需时间及速度增量
        microthrust: 微推推力，单位\muN
        sate_mass: 卫星质量，单位kg
        '''
        self.__microthrust = microthrust * 1.e-6
        self.__sate_mass = sate_mass

        self.__max_acceleration = self.__microthrust / self.__sate_mass

    def linear_motion(self, start_point, end_point):
        '''
        L = 0.5 * a * t^2
        "Bang-bang"机动模式
        '''
        distance = np.linalg.norm(end_point-start_point)
        delta_t = np.sqrt(self.__max_acceleration * distance)
        delta_v = delta_t * self.__max_acceleration

        return delta_t, delta_v

    def circling_motion(self, start_point, period=np.pi*2., repetition=1):
        from scipy.spatial.transform import Rotation
        rotate = Rotation.from_rotvec(period * np.array([0., 0., 1.], dtype=np.float64))
        end_point = rotate.apply(start_point)
        delta_t, delta_v = self.linear_motion(start_point, end_point)

        if repetition != 1:
            rotate = Rotation.from_rotvec(period * repetition * np.array([0., 0., 1.], dtype=np.float64))
            end_point = rotate.apply(start_point)

        return delta_t*repetition, delta_v*repetition, end_point

def get_max_JacobiSN():
    sn = mp.ellipfun('sn')
    def sn_scalar(u,m):
        return -np.abs(sn(np.linalg.norm(u), m).real)

    from scipy.optimize import minimize
    #sn_result = minimize(fun=sn_scalar, x0=(1.31), args=(-1), method = "Nelder-Mead", bounds=((0., 2.4),))
    sn_result = minimize(fun=sn_scalar, x0=(1.31), args=(-1), method = "L-BFGS-B", bounds=((0., 2.4),))
    return sn_result.x[0]

class circling_motion(object):
    def __init__(self, microthrust, sate_mass, max_angular_v, period=np.pi*2.) -> None:
        self.__microthrust = microthrust * 1.e-6
        self.__sate_mass = sate_mass
        self.__jacobi_max_val = get_max_JacobiSN()
        self.__max_angular_velocity = max_angular_v # arcsec/s

    def initialize(self, radius, period=np.pi/2.):
        self.__radius = radius
        self.__period = period
        self.__jacobi_parm = np.sqrt(np.power((self.__microthrust / (self.__sate_mass * self.__radius), 2.)))

    def angular_velocity(self, delta_t):
        def jacobi_integral(delta_t):
            return mp.ellipfun('sn', self.__jacobi_parm*delta_t + 0.j, -1)
        from scipy.optimize import fsolve
        pass





def main():
    #clr_clz = circling_motion(100., 1000., 100.)
    #from scipy import integrate
    #v, err = integrate.quad(clr_clz.angular_velocity, 0., 100.)
    #v = clr_clz.angular_velocity(100.)
    #print(v)
    print(get_max_JacobiSN())

if __name__ == "__main__":
    main()